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We investigate the effects of heterogeneous delays in the coupling of two excitable neural systems. 
Depending upon the coupling strengths and the time delays in the mutual and self-coupling, 
the compound system exhibits different types of synchronized oscillations of variable period. 
We analyze this synchronization based on the interplay of the different time delays and support 
the numerical results by analytical findings. In addition, we elaborate on bursting-like dynamics 
with two competing timescales on the basis of the autocorrelation function. 
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1. Introduction 

A plethora of synchronization phenomena has been found for coupled nonlinear oscillators in physical 
chemical and biologi cal systems |Pikovsky et al. , 2001; Boccaletti et al. , 2002; Mosekilde et al., 2002 



Balanov et al., 2009; Zhang et al., 201 1|. Especially the interplay of synchronization and time delay in 



coupled systems has received much interest recently [Scholl fc Schuster , 2008; Just et al., 2010; Atay 



2010 . Delayed coupling plays also a crucial role in the case of oscillation death |Choe et aT , 2007; Zou &; 



Zhan, 2009; Zou et al, 2012 1 and adaptive control schemes with delayed feedback applied to both chaotic 



[Wang et a/. j2010bl|Lehnert et a/.l|2011b| and non-chaotic system s |Selivanov et a/.[|2012|. Previous s tudies 



involved networks of a large number of elements |Atay et al. , 2004; Dhamala et al. , 2004; Kinzel et al. , 2009 
Choe et al.j |2010l |Zigzag et al.\ [20091 |Englert et al.[ |2010[ |Batista et al.j |2010HKanter et a/4|2011t|Flunkert 



&; Scholl, 2012 as well as simple recurring substructures consisting of a few systems only, so-called network 
motifs [Hauschildt et al\ [20061 P^ e et al \ [2007t |D'Huys eTaL " 



2008; Hovel et al, 2010a; Fiedler et al. 
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2010 



et al 



Flunkert et al , 2009 



2011; Adhikari et al 



Bra ndstetter et al\ |20T0l |D ? Huys et al] [20TT| |Hicke et al\ [20TT| |Kyrychko| 



2011 



. Considering the dynamics on networks with delay, the local elements can 



be either time-continuous or time-discrete as for iterated maps [Wang et al!| |2008a[ |2009b[. In addition 
a number of universal model-independent results have been obtained |Flunkert et al, 2010; Heiligenthal 



eT^|20lI] 



Synchronous dynamical patterns play also an important role in neuroscience 
Wang fc Lu[ |20051 |Hauptmann et al\ [20071 |Masoller et al 



et al 



2009; Senthilkumar et al, 2009; Liang et al, 2009; Lehnert et al 



20081 Wang et al 



2011a 



2008b, 2009a; Masoller 



Rossoni et al, 2005 



Popovych et al, 2011 



where on the one hand they can be observed in ensembles of neurons as pathological states like migraine, 
Parkinson's disease, or epilepsy. On the other hand synchronization can also be beneficial for recognition, 
learning, or neural information processing. Obviously the signal transmission between neurons in different 
brain areas is not instantaneous. Thus non-zero transmission times have to be taken into account as crucial 
quantities that influence the dynamics of individual neurons to a large extent. Consequently effects due 



to time delays have attracted more and more attention in the studies of neural networks Rossoni et al 



2009, 2010b 



Lehnert et al, 2011a; Kanter et al, 2011 Wang et al, 2011b and particularly in motifs of two couplec 



20051 IHauptmann et at] [2007| [M asolle r et al\ [20081 [FYiedrich fc Kinzelj [2009} |Wang et al4 |2010a[ |20lTa 



neuro ns |Scholl et al , 2009; Dahlem et al, 2009; Hovel et al, 2010a; Brandstetter et al, 2010; Hovel et al 



The latter can be seen as the smallest entity in a larger network. Interestingly, phenomena 
observed in this area of research show a strong similarity with recent findings in optoelectronic oscillators 
[Rosin et aZTj |2011| . 

Most previous works have assumed equal delay times in all connections. The focus of this paper, 
however, is on heterogeneous delays, which introduce additional timescales to the compound system. For 



this we consider a simple example of a network motif |Hauschildt et al , 2006; Dahlem et al , 2009; Panchuk 



et al, 2009; Hovel et al, 2010a|, i.e., two delay-coupled neurons with delayed self- feedback, and we assume 



all delay times to be different. This configuration might as well be understood as two effective populations 



of larger clusters of neurons with delayed internal and mutual connections |Vicente et al , 2008] . 

The rest of this paper is organized as follows: Section [2] introduces the neural model and the delay- 
coupling configurations. In Sec. [SJ we study interaction involving identical self-coupling delays numerically 
and analytically. The results are extended to the case of nonidentical self-coupling delays in Sec.[4j Section[5] 
considers bursting dynamics. Finally, we close with a conclusion in Sec. [6| 



2. Model 



We study a compound system of two coupled neural elements each represented by a FitzHugh-Nagumo 
model |FitzHugh , 1961 ; Nagumo et a/.[[l962] . The respective dynamic equations are paradigmatic for neural 
systems of type-II excitability, when periodic oscillations are generated in a Hopf bifurcation. We consider 
the case where the two elements are coupled such that each neural oscillator is subject to the delayed 
response from the other one. See Fig. [I] for a schematic diagram, where the time delays are denoted by 



r-C 



and 



r-C 



1 2 and C is the coupling strength. In addition, we take also delayed self- feedback Pyragas 
with a delay time r/^, i = 1, 2 and feedback strength K into account. 



1992 




Fig. 1. Schematic diagram of two coupled neural elements including the parameters of the mutual coupling (time delays 
T i i T 2 an( l coupling strength C) and the self-feedback (time delays and feedback gain K). 



The two coupled FitzHugh-Nagumo systems are described by the following set of delay differential 
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equations: 



ei x x = xi - y - yi + C [x 2 (t - t%) - xi(t)] 

+K[x 1 (t-rf c )-x 1 (t)] (la) 
yi=x 1 + a (lb) 



S 2 X 2 = x 2 - y - y 2 + C [xi (t - rf) - x 2 {t)] 

+K[x 2 {t-r^)-x 2 {t)] (lc) 
y 2 = x 2 + a, (Id) 

where e% denotes the timescale ratio between the slow inhibitor variable y\ and the fast activator variable 
Xi (i = 1, 2). The parameter a is known as threshold parameter. For \a\ > 1, the uncoupled system operates 
in the excitable regime. 



As it was shown in Ref. Panchu k et al\ [20091, two coupled systems with asymmetric delay times r-p, 
Tp can be reduced to a system with symmetric delay times r c . The difference between r-p and Tp leads 
only to a phase shift between oscillator 1 and 2. Therefore, we assume them to be equal rf = r§ = r c 
without loss of generality and Eqs. ([!]) can be rewritten as follows 

x 3 

e±i = x x - y - yi + C [x 2 (t - r c ) - xi(t)] 

+K [x! (t - t?) - X!(t)] (2a) 
y x = xi + a (2b) 

x 3 

ex 2 = x 2 - y - y 2 + C [xi (t - r c ) - x 2 {t)] 

+K[x 2 {t-T^)-x 2 {t)] (2c) 
V2 = x 2 + a. (2d) 

Throughout this paper, we choose the following set of parameters, unless specified otherwise: e\ — e 2 = 

£ = 0.01, a = 1.3, r c = 3, and C = 0.5. 

For \a\ > 1, the fixed point is always linearly stable (cf. Ref. |Panchuk et al |2009| ) and the system 



Eq. ([2]) shows various regular spiking and bursting patterns. Furthermore, several stable solutions can 
coexist for the same parameter values entailing high-level multi-stability. 

Before exploring the interplay between three different delay times, namely the mutual coupling delay r c 
and two nonidentical self-coupling delays ^ Tvf" (Sees, [i] and [5J, we will consider the case t± = = r K 
in the next Section. Analytical conditions for coherent spiking will be derived and generalized in the 
subsequent sections. 

3. Identical self-feedback delays 

The dynamics of the compound system ^ is diverse and hence, we will first introduce a classification. For 
this purpose, we will use an approach based on the mean interspike interval (ISI) \T^) = T^ l \ where 

N 



fjiWl j s the set of N interspike intervals for the time series of the z-th neuron (i = 1, 2) Dahlem et al. 

I J J 7=1 I i 

[20091 IScholl ~etll\ [2009] . A measure based on the ISI is a powerful tool to characterize regular, coherent 



spiking that is similar to a period- 1 orbit. Therefore, we take the ISI values for those cases into account 
that have small standard deviation. Failure of such a measure, e.g., for bursting dynamics, will be discussed 
in later sections. 

Figure [2] shows exemplary time series of coherent spiking for different values of K and r K , while the 
mutual coupling delay and strengths are fixed at r c = 3 and C = 0.5, respectively. The parameters are 
chosen as (K = 0.05, r K = 3), (K = 0.5, r K = 3), (K = 0.5, r K = 2), and (K = 0.5, r K = 4) in panels (a)- 
(d), respectively. All these combinations of r K and K exhibit coherent spiking, where the ISI is constant, 
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Fig. 2. Time series of Eqs. |2j) for identical self- feedback delays — = . The activator Xi and inhibitor yi are shown 
by blue (full) and red (dashed) curves, respectively. The parameters are chosen as (a) K = 0.05, r K = 3, (b) K = 0.5, r K = 3, 
(c) K — 0.5, r K — 2, and (d) K — 0.5, r K — 4. The black and light blue arrows show excitations due to self-feedback and 
mutual feedback, respectively. Other parameters: s = 0.01, a — 1.3, r C — 3, and C — 0.5. 



hence the standard deviation vanishes. Since the fixed point in the individual subsystems is stable in the 
excitable regime, we choose initial conditions such that only one subsystem is located in the fixed point 
whereas the other is subjected to a one-time-only excitation. This initial excitation eventually remains in 
the compound system due to the delayed coupling. 

Comparing Figs. [2^ a) and (b), one can see that small self- feedback gains, e.g., K — 0.05 in Fig. [21(a), 
are not able to trigger superthreshold excitations. Only subthreshold oscillations occur after times r as 
indicated by a black arrow in Figs. [2](a). In the following, we will derive analytical conditions for the delay 
times r K and r c such that regular, superthreshold spiking occurs. 

If the coupling strengths K and C are large enough, a spike at time t in one system will induce spikes 
at times t + r c in the other system (see light blue arrows) and t + r K in the first system (black arrows) by 
mutual coupling and self-feedback, respectively. The spike in the second system returns after a round-trip 
time 2t c . Thus, we have excitation events at times t + r K and t + 2r c in each subsystem. The spikes 
induced by these two sources of excitation become coherent, i.e., in resonance, when the delay time due to 
a round trip to the other system and back again, i.e., 2r c matches with the self-coupling delay r K . The 
same argument holds for integer multiples of r K and 2r^, leading to the following condition 

N K r K = N c 2t c , (3) 

with integer numbers N K and N c . 

Using this notation, Fig. [2] displays some combinations of N K and N c for coherently spiking states. 
Panels (a) and (b) refer to time delays tk = 3 = rc with N c = 1 and N K = 2. While panel (a) shows 
subthreshold oscillations after r K and superthreshold oscillations only after 2t c , panel (b) corresponds to 
a resonance yielding only fully pronounced, regular oscillations. Other values of the self-feedback delay, 
e.g., tk = 2 and tk = 4, result in different combinations of N K and N c . See, for instance, panels (c) 
and (d) that correspond to N K = 3, N c = 1, and N K = 3, N c = 2, respectively. For superthreshold 
oscillations, Eq. ^ yields the following condition for coherent spiking: 

N K ' W 



Synchronization of coupled neural oscillators with heterogeneous delays 5 



The corresponding ISIs are given by 



T 



2r 



c 



N K 



(5) 



on r 



with minimal integer numbers N K , TV^, i.e., the fraction N K /N c is irreducible. 

Equation ^ reflects the resonance condition N K spikes are induced by self-coupling during the 
round trip time 2r c of the mutual coupling. See Figs. j2^b)-(d). Figure [3] compares this analytical result to 
the numerical simulation. The green bars (positive T) show the numerically simulated ISI in dependence 
K for standard deviations smaller than 0.01. The other parameters are fixed at K = 0.5, C — 0.5, and 
3. For larger standard deviations the spiking is not coherent anymore. The red bars (negative T, i.e., 
inverted to facilitate comparison) show the analytically calculated delay times r K for which spiking occurs, 
and the value of the associated ISI using Eqs. Q, ^ and Q, (10). The latter describe the width of the 
resonance bars and will be derived later. It can be seen that the analytic results are in good agreement 
with the time delay t k , for which coherent spiking is found with the corresponding ISI T in the numerical 
simulation. Due to finite numerical accuracy, not all analytically possible ISI are detected in the time series. 
For large integers N c and N K , coherent spiking does not occur. 




Fig. 3. Period of possible ISIs of coherent spiking in dependence on tk for K — 0.5 and C — 0.5. The mutual time del ay i s 
fixed at r C = 3. Green: numerical simulation; red inverted values: analytical calculation using Eqs. ^ and Eqs. ([9]), |lo| ). 
The parentheses refer to some exemplary values (N C , N K ). Other parameters as in Fig. [2] 



It is possible to derive a condition concerning a possible phase shift of spikes in the activator variables 
x\ and X2 in the regime of coherent spiking as displayed, for instance, by Figs. [2^c) and (d). If we find 
integers N K and N c with 

N K r K = N c r c , (6) 

spikes in the first and second oscillator coincide. Thus for N K — N K /2 E N leading to N K even, there 
is no phase difference, i.e., we observe in-phase oscillations. Otherwise, for odd N K , the phase shift is tt 
corresponding to anti-phase oscillations. Using Eq. ^ we are able to predict for which values of r K and 
r c anti-phase oscillations occur. 

In order to introduce some helpful notation for the following derivation, Figs. Qa) andQb) display a 
time series and the respective phase portrait of a typical behavior of the neural system under consideration. 
A full excursion in phase space consists of a round trip from A through points B, C, and D back to A. 
The times for the transitions from A to B and from C to D are negligible, since the activator variables x\ 
change much faster than the inhibitors yi, due to timescale separation e <C 1. The transition from B to 
C happens close to the right slow branch of the cubic nullcline during a fixed firing time Tf = Tb^c- I n 



an earlier publication Scholl et al , 20091 we derived the following analytical approximation for the firing 
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with x\ 
ofT,* 



a + l 2 

= 2 and #i = 1 as an approximation for points B and C, respectively. Thus, we obtain a value 
0.45 for a = 1.3. Note that this approximation is valid for parameter values of a close to 1. For 
an improved estimate of Tf further away from the bifurcation point, one can use the position of the fixed 
point, i.e., intersection of the nullclines. For details see Appendix [I] For the final transition from D to A 
there remains the time Tj^^a — T — Tf, which will be considered later in this Section. 

In Fig. [3] one sees that each area of coherent spiking has a certain width. In the following we will derive 
an expression for the width Ar K of these coherence tongues by using the quantity Tf. For this, we will 
soften the condition ^ to within a certain tolerance: If the time shift \N K r K — N c 2r c \ is smaller than 
half the firing time T//2, spiking still occurs even though Eq. ^ is only approximately fulfilled. This leads 
to the following relation 



N 



K 



T K ± 



At k 



N c 2t c 



< 



(8) 



Here At k , i.e., the width of the coherence tongues, acts as tolerance in the timing of the incoming exci- 
tations. Equation ^ yields an upper bound for At k 

T f + 2\N C 2T C 



A- 



-K\ 



< 



N K r K \ 



N K 



Recalling condition Q, this simplifies to 



I At 



K\ 



< 



±1_ 

N K ' 



(9) 



(10) 



For a better analysis of the behavior of the compound system, it is helpful to investigate the ISI. 
Figure [5] shows the ISI as color code in the (if, r x )-plane for fixed r c — 3. Note that only spiking with 
an ISI standard deviation smaller than 0.01 is depicted. In the white region the standard deviation of the 
ISIs is large. Thus, white color marks the region, where no coherent spiking occurs. The horizontal lines 
correspond to delays r K in the self-coupling which are in resonance to r c . Some combinations of K and 



r are marked by black dots (a) - (d) referring to the time series in Fig. [2j 

The bright yellow region at small K and r K values refers to a regime, where oscillations with T = 
2t c = 6 dominate the dynamics due to mutual coupling, while the self-coupling leads only to subthreshold 
oscillations. In these cases, the self-coupling is too weak or too fast, i.e., small K or small r x , respectively, to 
initiate additional spikes. Compare also Fig. [2^a). The black region for large K and small r K corresponds 
to oscillation death due to the refractory phase of the neural oscillator |Scholl et al t 2009 . There the 
subsystem is not susceptible to an incoming activating signal. 
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Fig. 5. Interspike intervals of the spiking state obtained from numerical simulations in dependence on K and r with fixed 
r C = 3 and C = 0.5. The values marked by (a) to (d) refer to time series displayed in Fig. [2] Other parameters as in Fig. [2] 
A spike is considered as an excursion in phase space with x\ > 0. 



One can also calculate analytically the threshold at which the coherent spiking with T — 2r c ceases, 
i.e., the border between the bright yellow and the white regime in Fig. [5j At this boundary, self-coupling 
becomes strong enough to excite superthreshold spikes. To calculate an analogous boundary in Fig. [6j we 
set K — and vary C and r c and extend the analytic result later to nonzero self-coupling strength K to 
calculate the boundary in Fig. [5j Using the notation introduced above, the transition from D to A in Fig. [I] 
completes a full excursion. This last transition happens during a time interval Tj^^a — T — Tf. At this 
stage of the spike, the neural system is susceptible for the next excitation. For long periods T, the system 
relaxes to the fixed point xpp. For periods T used in this paper, however, the next excitation happens 
already at an earlier point A = (xi 5eru i, yi^nd)- An analytical estimate for Tf — Tb^c is given above, see 



Eq. ( 7b ) . We can derive a similar formula for T^^a 



r x l,end 1 rp^ 

T D ^ A = / : ~dxi . (11) 



If the system exhibits spikes with a period T = 2r c , we have T^^a — 2r c — Tf. Using this value, Eq. ( 11 ) 
yields an implicit expression for xi^ en d 

fp-2r c + T f \ , N 

Xl,end = - 2) exp I a 2 _i ) a ( 12 ) 

with the abbreviation p = a(2 + x^ end ) + 2-x\ end /2 that depends upon X\ en d- Since the relaxation from 
D to A follows closely the cubic y-nullcline, we can calculate a value for yi, en d as follows: 

3 

l,end / 10 \ 

Vl,end = Xl,end ' 



The mutual coupling (K = 0, C ^ 0) serves as an input in Eqs. pa] ) and ( |2c| ). Thus, it leads to a 
temporary, vertical shift of the cubic ?/-nullcline. If the minimum {xi imini yi, m i n ) of this dynamic ^/-nullcline, 
which includes input from the coupled system, is shifted beyond the point A= end? yi,end)i a spike is 
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triggered. The equation for the dynamic y-nullcline (for K 

3 

yi = ^i + y + C[(x 2 (t- 

r-C 



0) is given by 



(14) 



which depends on the delayed activator variable x 2 {t 
minimum of this nullcline can easily be calculated as 



r°) that eventually induces the next spike. The 



X\ . 



J l,min 



+ C[x 2 (t-T C )-X^ 



(15a) 
(15b) 



Finally, the condition y\^ en d = yi,min determines the bou ndary between coherent spiking and the quiescent 
state. Note that the delayed response x 2 (t — r c ) in Eq. (15b) remains to be chosen. 

The analytically calculated excitation threshold is shown in Fig. [6] as a dashed green curve for x 2 {t — 
r c ) = 1. This value is motivated by the assumption that (x 2 ,y 2 ) is located at point C in Fig.^b) at time 
t — r c . The ISI is depicted in color code. The black region refers to the quiescent state. 

Following the derivation described above, one can also calculate the excitation threshold for K ^ 
and fixed mutual coupling parameters C, r c in a similar way. For this, Eq. (14) needs to be extended as 
follows: 

x 3 

yi = Xl + ^ + C [(x 2 (t - r c ) - Xl ] 

-K[(x 2 (t-T K )- Xl ] , (16) 



Similarly Eqs. (15) become 



= -Vl-C-K 



+ K[x 2 (t 



J l,min 

3 



+ C [x 2 (t - T ) - X lymin ] 



(17a) 



(17b) 



Here one has to choose appropriate values for x 2 {i — r c ) and x 2 {i — r K ). 

Now we assume that the self-coupling delay r K sets the period of the regular spiking leading to 
Td^a = t k — Tf. Accordingly, the time 2r c in Eq. (12) of xi imin has to be replaced by r K . Then the 
condition derived from yi^ en d — 2/i, ? 



i.e., Eqs. (pi) and (fl7b|), and Eqs. (Tsl) and (flOl) for the ISI lead 
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7. Analytical calculation of the interspike intervals of the spiking state with r C = 3 and C = 0.5 using Eqs. Q, ^ and 



(|9), (10). The delayed values X2(t — r C ) and X2(t — r^) are chosen as —1.3 and 2. Other parameters as in Fig. [2] 



are 



to Fig. [7] that displays the analytically calculated ISIs. The delayed values x 2 {t — r c ) and x 2 (t — r 
chosen as —1.3 and 2, i.e., as points A and B in Fig. Qb), respectively. This figure is in good agreement 
with Fig. [5j which shows the numerically simulated ISIs. Not only the location of the coherence tongues 
are reproduced by the analytical formulas, but also their widths are in good agreement. 



4. Nonidentical self-coupling delays 

As one can see from the previous analysis, the dynamics of system ^ exhibits a variety of different solutions 
already for the simplifying restriction — = r K . In what follows we will consider 7^ r^, in which 
case one can expect even richer dynamics. The initial conditions are kept as before, i.e., one subsystem is 
initialized with a one-time excitation. 

In Fig. pi ISI diagrams in the (K, r/^) -parameter space are plotted for C — 0.5, r c — 3, and — 2 in 
color code. Again, white areas correspond to those solutions for which the ISI standard deviation is larger 
than a threshold value set to 0.02. The colored areas refer to parameters, for which coherent oscillations 
similar to a period- 1 orbit appear, associated with a single spike during one period. 

Comparing Fig. [8] with Fig. [5j i.e., the case of equal self- feedback delays, one finds common features 
like resonances or a boundary of the T — 2r c periodic dynamics at small K. A closer look, however, reveals 
also distinct differences. Some of the parameter regions of coherent spiking appear as tongues at specific 
ratios of the delays. This time, the regions with t± — and t± — 2r^ are more pronounced than for 
Ti — r c . The reason is that now the fixed self-coupling delay sets the basic timescale. 

Using the same argument as in Sec. [5] we can derive resonance conditions similar to Eq. Q. Namely, 
we need to require resonance for all pairs of the coupling delays: (i) t± and 2t c , (ii) and 2t c , as well 
as (hi) Ti and r 2 • The first two assumptions yield, as before, 

JVf 2r c = iVf rf (18) 

with irreducible integers iV-f /N± , and 

N$2t c = N*t? (19) 



with irreducible integers N% /Nrf • Furthermore, Eq. (18) divided by Eq. (19) immediately leads to 

7V lT f = N 2 7f (20) 

with Ni = N^N^/d and N 2 = N^N^/d, where d is the greatest common divisor of N^N^ and N^N^. 
Therefore, if the relations (18) and ( p~9| ) are satisfied, system ^ performs coherent spiking. Note that for 
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Fig. 8. Interspike intervals T in the yK, t± J -plane with r — 3,r 2 = 2 as color code. Only solutions with an ISI standard 
deviation smaller than 0.02 are shown. Other parameters as in Fig. [2] 



any three rational numbers t c , r/^, and r 2 there always exist corresponding integers iV-p , , TVp, N^f , 
Ni and N 2l for which the Eqs. @-([20) hold. 

The period of the coherent solution is also obtained by analogy with the previous case. First, we find 
from Eqs. (ph and dl9h 



Ti 



T 2 



2r 



c 



Finally using Eq. (20) we have 



Thus, the estimated period is 



N* 



± 



N\ N 2 



(21) 
(22) 

(23) 
(24) 



T = mm{T 1 ,T 2 ,T 3 }. 

As a consequence, if Nf, Nf, N ,^, and N 2 are large, the estimated oscillation period becomes small. 
Then the coherent solution cannot be realized due to the refractory phase of the neural subsystems. 

In order to determine the width of the areas for coherent spiking in the (K, t^*) -plane, we use a similar 
reasoning as before. From Eqs. (18) and (20) we get 



and 



Ni 



K 



At, 



K 



N c 2t c 



< 



N 2 T~2 



K 



< 



2 ' 



respectively. As a result, the width of the regular spiking regime simplifies to 

|Arf|<min{^,^f 



(25) 



(26) 



(27) 



Figure [9] shows the analytical estimate for the tongues of coherent spiking corresponding to Eqs. ( [18] ) 
and (19). The region for the solution of period T « 2r c is obtained numerically at small feedback gains 
K . There is a good correspondence to Fig. [51 however, some approximated regions for solutions with small 
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Fig. 9. Analytical approximation of tongues for coherent spiking with period less than 2t C in the (k, t^j -plane (t C = 

3, = 2). The firing time Tf is set to 0.38 (see Appendix^). The yellow region for the solution of period 2t G at small K is 
calculated numerically. Other parameters as in Fig. [2] 



periods are wider than the simulated ones. The reason is that the approximation for the firing time XV, 
given by Eq. (A.l), is still too rough. For small periods, the solution looks as exemplarily depicted in 



Fig. 10, and the derivation for TV does not hold anymore. 
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Fig. 10. Time series and phase portrait of the first subsystem variables x\ and y\ as blue and red curves for the periodic 
solution with period T w 0.5. Self-coupling parameters: K — 0.5, — 0.5, — 2. Other parameters as in Fig. [2] 



As a generalization of our previous analytical approach, we can derive a relation between all three delay 
times in general form. Provided that the parameters C and K are large enough to yield superthreshold 
excitation, the spike of the first neuron at time t will induce spikes at times t+r^ and t + 2r c . Similarly the 

second neuron will spike at time t + r c , as well as at times t + r^ + Tvf" , t-\-r c + , t + T C + 37vf" , Thus, 

the first subsystem will get also excitation impulses at times t + 2r c + Tvf" , t + 2r c + 27vf" , t + 2r c + 3t^", 
and so on. From this we get 

mirf = h(2r c + mrf ), (28) 

where l\, mi, 724 are arbitrary positive integers. The same discussion is applicable for the second subsystem, 
and therefore we can also write down the symmetric condition 

n 2 T? = / 2 (2T C + m 2 rf) (29) 
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with positive integers fa^m^^n^. Adding Eqs. (|28j) and ([29]) we derive 

2(/i + fa)r C = (mi - fam 2 )rf + (n 2 - l\n x )r^ . 

Denoting / = 2(/i + fa)/*!, m — (mi — fam^/d, and n = (712 — l\n\)/d, where d is the greatest common 
divisor of the three numbers 2(/i + fa), mi — farri2, and 712 — Zini, we get 



rfiTi + TIT2 



(30) 



with integers / > 0, and m, n being of any sign. Note that Eq. (30) can also be obtained directly from 



adding or subtracting Eqs. (18) and (19). 



At last, by moving all the terms onto one side of the equation and relabeling, Eq. (30) can be rewritten 
in the general form: 



(31) 



lr c + mrf + TIT2 = 0, 

where /, m, n are arbitrary integers of any sign. It should be mentioned that the relation similar to Eq. ( |3l| ) 
was obtained in Ref. Zigzag et al. (2009 for a two-dimensional time-discrete system with several non-equal 
delays. 




If we fix t c , Eq. (31 ) defines different lines in the (r-j^, r^)-plane, depending on /, m, and n. Figure 11 ^a) 
visualizes multiple combinations as black lines on top of the ISI T shown in color code. The regions, in 
which periodic solutions were found numerically, accumulate along the lines that are given by Eq. (31). In 
Fig. [TT|b) these analytically obtained conditions are separately shown including the values of /, m, and n. 
The analytical results are in excellent agreement with the numerical calculations. 



5. Bursting and autocorrelation function 



In Fig. [TT[a) of the previous section, the lines obtained analytically from the resonance condition (31) agree 
very well with the numerical results on the ISI of regular spiking. However, the regions of periodic spiking 
fill only a part of the parameter plane. The ISI approach does not allow for an analysis of bursting-type 
solutions. There, the time series exhibits bunched spiking patterns and thus different timescales. 

Therefore, we consider the autocorrelation function (ACF) ^(s) as an alternative tool for the analysis 
of the coupled system dynamics |Hauschildt et al. , 2006| . The ACF ^(s) of an arbitrary time series x{t) is 
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Fig. 12. Maxima of the autocorrelation function ^ in the ^rf^, -plane for 
given by the Eq. (30). Other parameters as in Fig. [5] 



3,K — 0.5. The black lines lines are 



defined as 



([x(t 



(x)][x(t) - (x)]) . 



(32) 



where the averages (•) are calculated over the whole simulated time interval. The value a denotes the 
standard deviation of x{t), i.e., a 2 = ([x(t) — (x)] 2 ). 

To determine whether a given solution shows quasi-periodic neural activity, we calculate the ACF *f?(s) 
of ^i(t), which has several maxima. The first, trivial maximum is, obviously, obtained at s = and equals 
1. The second largest maximum marks the best coincidence between the original and the shifted time 
series, and is found for 5 = 5*, which characterizes the length of the repeated (periodic or quasi-periodic) 
pattern, and in case of periodic spiking equals the ISI. 



The result of this alternative analysis based on the ACF is depicted in Fig. [12} In this figure, the 
resonances of Eq . (|3l| ) (Fig. [TTj^b)) are fully visible and well separated from each other, cf. the yellow and 
red areas in Fig. 12 that are missing in Fig. [TTla). 



The ACF allows for investigation of both periodic and quasi-periodic behavior. Figure [13] shows time 
series and phase space projections for bursting-type solutions. Since the time series consists of repeated 
bunches of spikes, the mean interspike interval is no longer a good measure. This is due to two different 
timescales in the activator variables x\ and X2- To investigate this bursting- type behavior we use again the 
ACF analysis. 

Figure [l4^b) displays the ACF of the time series from Fig. 13, In Fig. [l4^a) the time series for the 



first activator x\ is shown again for convenience and Fig. 14 [b) presents the ACF It can be seen 

that the ACF approaches unity for the second time - after the trivial perfect correlation at s = - at a 
displacement of s = s* « 2.01, which equals the repetition period of the bursting pattern. The bursts can 
also be resolved by the ACF as the fast oscillations. To conclude, the ACF enables to distinguish between 
inter-burst and intra-burst timescales. 



6. Conclusion 

We have investigated effects of heterogeneous time delays for mutual and self-coupling of a simple network 
motif that consists of two neural systems. This setup is realized by two delay-coupled elements of FitzHugh- 
Nagumo type, which is a paradigmatic model of neural interaction. The two subsystems operate in the 
excitable regime, and excitation occurs due to the incoming delayed signals via both mutual and self- 
coupling. 

At first, we have considered identical self-coupling delays and analyzed the regular periodic dynamics 



14 A. Panchuk et al 




2990 3000 -2 2 -1 1 

(d) t (e) x 2 (f) Yl 



Fig. 13. Time series (a, d) and phase space projections (b, c, e, f) of xuyi and X2,V2 for the bursting- type solution of the 
system |2j) with K — 0.5, t± =2.2, — 2. Other parameters as in Fig. V2\ 




on the basis of the mean interspike interval. For small feedback strengths, the system exhibits regular 
behavior with the period of about twice the mutual coupling delay. With increasing feedback strength, 
however, the self-feedback term becomes stronger and the system can perform more frequent spikes due 
to this additional source of excitation. This happens if the mutual coupling and self-feedback delays are 
in resonance. In the parameter plane of self-coupling delay and strength, the regions where the regular 
oscillations exist, resemble stripes emerging from the resonance with the mutual couping delays. We have 
provided an analytical formula for these resonance conditions and the period of the synchronized oscilla- 
tions. A comparison with numerical simulations shows excellent agreement. 

Next, we have focused on the case of non-identical self-coupling delays. We have considered two- 
dimensional projections of the parameter space and measured the regularity of the dynamics by the mean 
interspike interval as well as by the autocorrelation function. Similar to the case of identical self-coupling 
delays, we observe synchronized periodic dynamics, if the three delay times satisfy a resonance condition. 
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We have also derived a formula for the period of the synchronized regular dynamics, and compared the 
theoretical results with our numerical simulations. 

Finally, we have studied bursts of subsequent excitation spikes. To identify such solutions the measure 
of interspike intervals is not appropriate any more, and we have used the autocorrelation function as an 
alternative tool for the analysis of the dynamics. 

In conclusion, we have shown that heterogeneous time delays give rise to regular synchronization 
patterns of different periods, which depend upon resonance conditions of the involved time delays. 

We have restricted our investigations to local dynamics of type-II excitability related to a Hopf bifur- 
cation (FitzHugh-Nagumo model). It is also interesting to consider other models describing, for instance, 
type-I excitability involving a saddle-node bifurcation on an invariant cycle or physiologically oriented 
models such as Hodgkin-Huxley- or Morris-Lecar-like equations, but those studies are beyond the scope of 
the present paper. Another important direction for future research is to increase the number of elements. 
First approaches to study delayed coupling in large networks have already been reported (see references in 
Sec. [I]). In principle, the results of the presented work can also be applied to coupled systems of more than 
two elements. Then one has to carefully analyze the combinatorics of the various exciting self-feedback 
and cross-coupling pluses based on the given network topology. In the case of networks, asymmetries in 
the coupling strengths can also become important, e.g., if multiple subthreshold excitations accumulate. In 
addition, different coupling strengths play a crucial role in the presence of both excitatory and inhibitory 
connections in neural networks. In the presented study of two coupled elements, however, incoming pulses 
result in all-or-nothing events: if the signal is large enough, it will trigger a full-scale excitation. Otherwise, 
only a small, subthreshold excitation is possible. 
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Appendix A Approximation of the firing time 



The formula (7b) was derived only for a close to unity, i.e., close to the bifurcation point. For larger 
values of a, however, a similar estimate holds. In fact, the idea is to integrate in time along the right branch 
of the cubic nullcline from the point B to the point C (See Fig. [4j). If a is close to 1, the coordinates of 
these points can be approximated by £>(2, —2/3) and C(l,2/3). For large a this approximation is rather 
bad. 

To improve the formula one can do the following. The x-coordinate of point A cannot exceed —a, 
because the cycle cannot cross the fixed point P(— a, a 3 /3 — a). Therefore, the maximum coordinates for A 
are approximately (—a, a 3 /3 — a). Since the transition from the left branch to the right branch of the cubic 
nullcline happens almost instantaneously, the y-coordinate of B could also be approximated as a 3 /3 — a. 
This yields B = (a/2 + y/12 - 3a 2 /2, a 3 /3 - a). 

To calculate the coordinates of the point C we proceed as follows. Since the periodic trajectory in 
Fig. |4^b) appears to be symmetric with respect to the origin (0, 0), the y-coordinate of C can be obtained 
as a — a 3 /3, and hence C(a, a — a 3 /3). This leads to the following improved approximation for the firing 
time T f 
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T f = / l -^dx x 
Jx B xi+a 

= - a(x B - x c ) + (a 2 - 1) In 

,o N , 3a + Vl2 - 3a 2 a / , -\ 3 

= (a 2 -l)ln ^- -(a+ Vl2-3a2)+- 

This yields a value of T f « 0.38 for a = 1.3. 
References 

Adhikari, B. M., Prasad, A. & Dhamala, M. [2011] "Time-delay-induced phase-transition to synchrony in 

coupled bursting neurons," Chaos 21, 023116. 
Atay, F. M. (ed.) [2010] Complex Time-Delay Systems, Understanding Complex Systems (Springer, Berlin 

Heidelberg). 

Atay, F. M., Jost, J. & Wende, A. [2004] "Delays, connection topology, and synchronization of coupled 

chaotic maps," Phys. Rev. Lett. 92, 144101. 
Balanov, A. G., Janson, N. B., Postnov, D. E. & Sosnovtseva, O. V. [2009] Synchronization: From Simple 

to Complex (Springer, Berlin). 
Batista, C. A. S., Lopes, S. R., Viana, R. L. & Batista, A. M. [2010] "Delayed feedback control of bursting 

synchronization in a scale- free neuronal network," Neural Networks 23, 114-124. 
Boccaletti, S., Kurths, J., Osipov, G., Valladares, D. L. & Zhou, C. S. [2002] "The synchronization of 

chaotic systems," Phys. Rep. 366, 1-101. 
Brandstetter, S. A., Dahlem, M. A. & Scholl, E. [2010] "Interplay of time-delayed feedback control and 

temporally correlated noise in excitable systems," Phil. Trans. R. Soc. A 368, 391-421. 
Choe, C. U., Dahms, T., Hovel, P. & Scholl, E. [2010] "Controlling synchrony by delay coupling in networks: 

from in-phase to splay and cluster states," Phys. Rev. E 81, 025205(R). 
Choe, C. U., Flunkert, V., Hovel, P., Benner, H. & Scholl, E. [2007] "Conversion of stability in systems 

close to a Hopf bifurcation by time-delayed coupling," Phys. Rev. E 75, 046206. 
Dahlem, M. A., Hiller, G., Panchuk, A. & Scholl, E. [2009] "Dynamics of delay-coupled excitable neural 

systems," International Journal of Bifurcation and Chaos 19, 745-753. 
Dhamala, M., Jirsa, V. K. & Ding, M. [2004] "Enhancement of neural synchrony by time delay," Phys. 

Rev. Lett. 92, 074104. 

D'Huys, O., Fischer, I., Danckaert, J. & Vicente, R. [2011] "Role of delay for the symmetry in the dynamics 

of networks," Phys. Rev. E 83, 046223. 
D'Huys, O., Vicente, R., Erneux, T., Danckaert, J. & Fischer, I. [2008] "Synchronization properties of 

network motifs: Influence of coupling delay and symmetry," Chaos 18, 037116. 
Englert, A., Kinzel, W., Aviad, Y., Butkovski, M., Reidler, L, Zigzag, M., Kanter, I. & Rosenbluh, M. 

[2010] "Zero lag synchronization of chaotic systems with time delayed couplings," Phys. Rev. Lett. 

104, 114102. 

Fiedler, B., Flunkert, V., Hovel, P. & Scholl, E. [2010] "Delay stabilization of periodic orbits in coupled 

oscillator systems," Phil. Trans. R. Soc. A 368, 319-341. 
FitzHugh, R. [1961] "Impulses and physiological states in theoretical models of nerve membrane," Biophys. 

J. 1, 445-466. 

Flunkert, V., D'Huys, O., Danckaert, J., Fischer, I. & Scholl, E. [2009] "Bubbling in delay-coupled lasers," 
Phys. Rev. E 79, 065201 (R). 

Flunkert, V. & Scholl, E. [2012] "Chaos synchronization in networks of delay-coupled lasers: Role of the 
coupling phases," New. J. Phys. 14, 033039. 

Flunkert, V., Yanchuk, S., Dahms, T. & Scholl, E. [2010] "Synchronizing distant nodes: a universal classi- 
fication of networks," Phys. Rev. Lett. 105, 254101. 



(A.l) 
(A.2) 
(A.3) 



REFERENCES 17 



Friedrich, J. & Kinzel, W. [2009] "Dynamics of recurrent neural networks with delayed unreliable synapses: 

metastable clustering," J. Comput. Neurosci. 27, 65-80. 
Hauptmann, C., OmePchenko, O. E., Popovych, O. V., Maistrenko, Y. L. & Tass, P. A. [2007] "Control of 

spatially patterned synchrony with multisite delayed feedback," Phys. Rev. E 76, 066209. 
Hauschildt, B., Janson, N. B., Balanov, A. G. & Scholl, E. [2006] "Noise-induced cooperative dynamics 

and its control in coupled neuron models," Phys. Rev. E 74, 051906. 
Heiligenthal, S., Dahms, T., Yanchuk, S., Jungling, T., Flunkert, V., Kanter, I., Scholl, E. & Kinzel, W. 

[2011] "Strong and weak chaos in nonlinear networks with time-delayed couplings," Phys. Rev. Lett. 

107, 234102. 

Hicke, K., D'Huys, O., Flunkert, V., Scholl, E., Danckaert, J. & Fischer, I. [2011] "Mismatch and syn- 
chronization: Influence of asymmetries in systems of two delay-coupled lasers," Phys. Rev. E 83, 
056211. 

Hovel, P., Dahlem, M. A., Dahms, T., Hiller, G. & Scholl, E. [2009] "Time-delayed feedback control of 

delay-coupled neurosystems and lasers," Preprints of the Second IFAC meeting related to analysis 

and control of chaotic systems (CHAOS09) (World Scientific), ( arXiv:0912.3395[ ). 
Hovel, P., Dahlem, M. A. & Scholl, E. [2010a] "Control of synchronization in coupled neural systems by 

time-delayed feedback," International Journal of Bifurcation and Chaos 20, 813-815. 
Hovel, P., Shah, S. A., Dahlem, M. A. & Scholl, E. [2010b] "Feedback-dependent control of stochastic 

synchronization in coupled neural systems," From physics to control through an emergent view, eds. 

Fortuna, L., Fradkov, A. L. & Frasca, M. (World Scientific, Singapore), pp. 35-44. 
Just, W., Pelster, A., Schanz, M. & Scholl, E. [2010] "Delayed complex systems," Theme Issue of Phil. 

Trans. R. Soc. A 368, pp.301-513. 
Kanter, I., Kopelowitz, E., Vardi, R., Zigzag, M., Kinzel, W., Abeles, M. & Cohen, D. [2011] "Nonlocal 

mechanism for cluster synchronization in neural circuits," Europhys. Lett. 93, 66001. 
Kinzel, W., Englert, A., Reents, G., Zigzag, M. & Kanter, I. [2009] "Synchronization of networks of chaotic 

units with time-delayed couplings," Phys. Rev. E 79, 056207. 
Kyrychko, Y. N., Blyuss, K. B. & Scholl, E. [2011] "Amplitude death in systems of coupled oscillators with 

distributed-delay coupling," Eur. Phys. J. B 84, 307-315. 
Lehnert, J., Dahms, T., Hovel, P. & Scholl, E. [2011a] "Loss of synchronization in complex neural networks 

with delay," Europhys. Lett. 96, 60013. 
Lehnert, J., Hovel, P., Flunkert, V., Guzenko, P. Y., Fradkov, A. L. & Scholl, E. [2011b] "Adaptive tuning 

of feedback gain in time-delayed feedback control," Chaos 21, 043111. 
Liang, X., Tang, M., Dhamala, M. &; Liu, Z. [2009] "Phase synchronization of inhibitory bursting neurons 

induced by distributed time delays in chemical coupling," Phys. Rev. E 80, 066202. 
Masoller, C., Torrent, M. C. & Garcia-Ojalvo, J. [2008] "Interplay of subthreshold activity, time-delayed 

feedback, and noise on neuronal firing patterns," Phys. Rev. E 78, 041907. 
Masoller, C., Torrent, M. C. & Garcia-Ojalvo, J. [2009] "Dynamics of globally delay-coupled neurons dis- 
playing subthreshold oscillations," Philosophical Transactions of the Royal Society A: Mathematical, 

Physical and Engineering Sciences 367, 3255-3266. 
Mosekilde, E., Maistrenko, Y. L. & Postnov, D. [2002] Chaotic Synchronization: Applications to Living 

Systems (World Scientific, Singapore). 
Nagumo, J., Arimoto, S. & Yoshizawa., S. [1962] "An active pulse transmission line simulating nerve axon." 

Proc. IRE 50, 2061-2070. 

Panchuk, A., Dahlem, M. A. & Scholl, E. [2009] "Regular spiking in asymmetrically delay-coupled 
FitzHugh-Nagumo systems," Proc. NDES 09 , 177-17E |aDriv:0911.2071 vl. 

Pikovsky, A. S., Rosenblum, M. G. & Kurths, J. [2001] Synchronization, A Universal Concept in Nonlinear 
Sciences (Cambridge University Press, Cambridge). 

Popovych, O. V., Yanchuk, S. & Tass, P. A. [2011] "Delay- and coupling-induced firing patterns in oscil- 
latory neural loops," Phys. Rev. Lett. 107, 228102. 

Pyragas, K. [1992] "Continuous control of chaos by self-controlling feedback," Phys. Lett. A 170, 421. 

Rosin, D. P., Callan, K. E., Gauthier, D. J. & Scholl, E. [2011] "Pulse-train solutions and excitability in 
an optoelectronic oscillator," Europhys. Lett. 96, 34001. 



18 REFERENCES 

Rossoni, E., Chen, Y., Ding, M. & Feng, J. [2005] "Stability of synchronous oscillations in a system of 
Hodgkin-Huxley neurons with delayed diffusive and pulsed coupling," Phys. Rev. E 71, 061904. 

Schoil, E., Hiller, G., Hovel, P. & Dahlem, M. A. [2009] "Time-delayed feedback in neurosystems," Phil. 
Trans. R. Soc. A 367, 1079-1096. 

Scholl, E. & Schuster, H. G. (eds.) [2008] Handbook of Chaos Control (Wiley- VCH, Weinheim), second 
completely revised and enlarged edition. 

Selivanov, A. A., Lehnert, J., Dahms, T., Hovel, P., Fradkov, A. L. & Scholl, E. [2012] "Adaptive synchro- 
nization in delay-coupled networks of Stuart-Landau oscillators," Phys. Rev. E 85, 016201. 

Senthilkumar, D. V., Kurths, J. & Lakshmanan, M. [2009] "Inverse synchronizations in coupled time- 
delay systems with inhibitory coupling," Chaos: An Interdisciplinary Journal of Nonlinear Science 
19,023107. 

Vicente, R., Gollo, L. L., Mirasso, C. R., Fischer, I. & Gordon, P. [2008] "Dynamical relaying can yield 

zero time lag neuronal synchrony despite long conduction delays," Proc. Natl. Acad. Sci. 105, 17157. 
Wang, Q., Chen, G. & Perc, M. [2011a] "Synchronous bursts on scale- free neuronal networks with attractive 

and repulsive coupling," PLoS ONE 6, el5851. 
Wang, Q., Duan, Z., Perc, M. & Chen, G. [2008a] "Synchronization transitions on small- world neuronal 

networks: Effects of information transmission delay and rewiring probability," EPL 83, 50008. 
Wang, Q., Lu, Q. & Chen, G. [2008b] "Synchronization transition induced by synaptic delay in coupled 

fast-spiking neurons," International Journal of Bifurcation and Chaos 18, 1189. 
Wang, Q., Lu, Q., Chen, G., Feng, Z. & Duan, L. [2009a] "Bifurcation and synchronization of synaptically 

coupled FHN models with time delay," Chaos, Solitons and Fractals 39, 918-925. 
Wang, Q., Perc, M., Duan, Z. & Chen, G. [2009b] "Synchronization transitions on scale- free neuronal 

networks due to finite information transmission delays," Phys. Rev. E 80, 026206. 
Wang, Q., Perc, M., Duan, Z. & Chen, G. [2010a] "Impact of delays and rewiring on the dynamics of 

small-world neuronal networks with two types of coupling," Physica A 389, 3299-3306. 
Wang, Q.-Y., Aleksandra, M., Matjaz, P. & Lu, Q.-S. [2011b] "Taming desynchronized bursting with delays 

in the macaque cortical network," Chinese Physics B 20, 040504. 
Wang, Q.-Y. & Lu, Q.-S. [2005] "Time delay-enhanced synchronization and regularization in two coupled 

chaotic neurons," Chinese Physics Letters 22, 543. 
Wang, Q. Y., Lu, Q. S. & Duan, Z. S. [2010b] "Adaptive lag synchronization in coupled chaotic systems 

with unidirectional delay feedback," Int. J. Nonlin. Mech. 45, 640-646. 
Zhang, W., Tang, Y., Fang, J. -a. & Zhu, W. [2011] "Exponential cluster synchronization of impulsive 

delayed genetic oscillators with external disturbances," Chaos 21, 043137. 
Zigzag, M., Butkovski, M., Englert, A., Kinzel, W. & Kanter, I. [2009] "Zero-lag synchronization of chaotic 

units with time-delayed couplings," Europhys. Lett. 85, 60005. 
Zou, W., Tang, Y., Li, L. & Kurths, J. [2012] "Oscillation death in asymmetrically delay-coupled oscilla- 
tors," Phys. Rev. E 85, 046206. 
Zou, W. & Zhan, M. [2009] "Partial time-delay coupling enlarges death island of coupled oscillators," Phys. 

Rev. E 80, 065204. 



